Thalamocortical Connections Mature along a Sensorimotor-Association Axis of Cortical Developmental Heterochronicty
Read in Data for Developmental Characterization
Glasser parcel names and mesulam assignments
#Glasser regions with corresponding labels, tract names, and mesulam assignments
glasser.regions <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360_regionlist.csv") #parcel name, label name
glasser.regions$tract <- paste0("thalamus_", glasser.regions$orig_parcelname) %>% gsub("_ROI", "_autotrack", .) %>% gsub("-", "_", .) #create tract variable with format thalamus_$hemi_$region_autotrack, no dashes -
#Glasser regions assigned to mesulam zones
glasser.assignments <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser360-mesulam_economo_yeo-assignments.csv")
glasser.assignments <- merge(glasser.assignments, glasser.regions, by = "label", sort = F)
glasser.assignments <- glasser.assignments %>% select(label, tract, orig_parcelname, mesulam_assignment)S-A axis
S.A.axis.cifti <- read_cifti("/cbica/projects/thalamocortical_development//Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii") #S-A_ArchetypalAxis repo
S.A.axis <- data.frame(SA.axis = rank(S.A.axis.cifti$data), orig_parcelname = names(S.A.axis.cifti$Parcel))
S.A.axis <- merge(S.A.axis, glasser.assignments, by = c("orig_parcelname"), sort = F)
S.A.axis$SA.axis.bin <- as.numeric(cut2(S.A.axis$SA.axis, g = 5)) #divide the S-A axis into 5 ranked bins, 72 regions per bin
S.A.axis %>% group_by(SA.axis.bin) %>% do(num.regions = length(.$SA.axis), mean.rank = mean(.$SA.axis)) %>% unnest(cols=c("num.regions", "mean.rank"))## # A tibble: 5 × 3
## SA.axis.bin num.regions mean.rank
## <dbl> <int> <dbl>
## 1 1 72 36.5
## 2 2 72 108.
## 3 3 72 180.
## 4 4 72 252.
## 5 5 72 324.
NeuroSynth cognitive atlas term maps
#Neurosynth z-score maps generated by nimare. Nimare uses multilevel kernel density analysis- Chi-square analysis + FDR-correction to use the same procedure as Neurosynth
neurosynth.terms <- read.csv("/cbica/projects/thalamocortical_development/Maps/neurosynth/atl-glasser360_neurosynth.csv") %>% select(-regionName) %>% select(-timing) #neurosynth term meta-analytic scores for 123 terms present in both NeuroSynth and the Cognitive Atlas, ordered in lh --> rh glasser order
neurosynth.termlist <- names(neurosynth.terms) #list of terms (cue willy wonka line "dont you want to know our names?" "cant imagine why it wouldnt matter")
neurosynth.terms$label[1:180] <- glasser.regions$label[181:360] #lh first
neurosynth.terms$label[181:360] <- glasser.regions$label[1:180] #rhSpatial permutation (spin) test parcel rotation matrix
perm.id.full <- readRDS("/cbica/projects/thalamocortical_development/software/rotate_parcellation/glasser_sphericalrotations_N10000.rds") #10,000 spatial null spins
spin.parcels <- glasser.regions #order of complete set of glasser parcels for spinningDataset-specific diffusion measure spreadsheets
FA.glasser.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_FA_finalsample_primary.csv") #generated by sample_construction/PNC/tractmeasures_dfs_PNC.R
FA.glasser.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_FA_finalsample_primary.csv") #generated by /sample_construction/HCPD/tractmeasures_dfs_HCPD.RDataset-specific tract lists
tractlist.pnc <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/PNC/PNC_tractlist_primary.csv") #generated by tract_measures/tractlists/thalamocortical_tractlists.R
tractlist.hcpd <- read.csv("/cbica/projects/thalamocortical_development/thalamocortical_structuralconnectivity/individual/HCPD/HCPD_tractlist_primary.csv") #generated by tract_measures/tractlists/thalamocortical_tractlists.RDevelopmental statistics
setwd("/cbica/projects/thalamocortical_development/thalamocortical_results/development_results/") #gam output dir
#list all FA development files in development results directory for accessing. Files were generated by gam_functions/fit_ageGams.R
files <- list.files(getwd(), pattern='FA', ignore.case=T, full.names = F)
#files <- files[!str_detect(files,pattern="temporalSAcorr")] #except these, which we will read in later because they are very large
filenames <- c()
#generate variable names to assign files data to
for(name in files){
Rname <- gsub('^.{12}|.{4}$', '', name) #remove ".csv" from end of filename and "development_" from start of filename
filenames <- append(filenames, Rname) #save filenames into a character vector
}
#read in files and assign to variables
for(i in 1:length(filenames)){
Rfilename <- sprintf("%s",filenames[i]) #save filename as an individual string
if(grepl("gameffects", Rfilename)) {
x <- read.csv(files[i], header = TRUE)
x <- merge(x, glasser.assignments, by = "tract")
assign(Rfilename, x)
}
if(!grepl("temporal", Rfilename) & !grepl("gameffects", Rfilename)) {
x <- read.csv(files[i], header = TRUE)
x <- merge(x, S.A.axis, by = "tract")
assign(Rfilename, x)
}
if(grepl("temporal", Rfilename)){
x <- read.csv(files[i], header = FALSE)
assign(Rfilename, x)
}
}
rm(x)A Spectrum of Thalamocortical Developmental Trajectories
Cortex-wide thalamic connectivity developmental profiles
Number of significant developmental effects
#Function to calculate the number and percent of thalamocortical connections showing a significant developmental change in a given measure
sig.effects <- function(measure, atlas, dataset){
ageeffects.df <- get(sprintf("gameffects_%s_%s_%s", measure, atlas, dataset))
ageeffects.df$significant <- p.adjust(ageeffects.df$Anova.smooth.pvalue, method = c("fdr")) < 0.05 #fdr-corrected significant connections
sigeffect.totaln <- sum(ageeffects.df$significant) #total number of significant connections
sigeffect.percent <- round(sigeffect.totaln/length(ageeffects.df$significant)*100, 2) #percent of significant connections
print(sprintf("In %s, %s thalamocortical connections (%s percent) show significant age-related changes in %s", dataset, sigeffect.totaln, sigeffect.percent, measure))
}PNC
## [1] "In pnc, 201 thalamocortical connections (90.13 percent) show significant age-related changes in FA"
HCPD
## [1] "In hcpd, 174 thalamocortical connections (75.65 percent) show significant age-related changes in FA"
Thalamocortical connection GAM smooth functions
PNC
smoothcentered_FA_glasser_pnc.plot <- merge((smoothcentered_FA_glasser_pnc %>% filter(grepl("thalamus_R", tract))), (gameffects_FA_glasser_pnc %>% select(tract, GAM.smooth.partialR2)), by = "tract") #zero-centered developmental smooth functions for RH connections + tract-specific partial R2 for plotting
ggplot(smoothcentered_FA_glasser_pnc.plot, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) +
geom_line(linewidth = .3, alpha = .8) +
scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(0, 0.1), oob = squish) +
theme_minimal() +
labs(x="\nAge", y="Connection FA\n") +
scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
theme(
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Developmental_Trajectories_PNC.pdf", device = "pdf", dpi = 500, width = 2.02, height = 1.65)HCPD
smoothcentered_FA_glasser_hcpd.plot <- merge((smoothcentered_FA_glasser_hcpd %>% filter(grepl("thalamus_R", tract))), (gameffects_FA_glasser_hcpd %>% select(tract, GAM.smooth.partialR2)), by = "tract")
ggplot(smoothcentered_FA_glasser_hcpd.plot, aes(x = age, y = est, group = tract, color = GAM.smooth.partialR2)) +
geom_line(linewidth = .3, alpha = .8) +
scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.02, 0.075), oob = squish) +
theme_minimal() +
labs(x="\nAge", y="Connection FA\n") +
scale_x_continuous(breaks=c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0,.45)) +
scale_y_continuous(limits = c(-0.02, 0.0138)) +
theme(
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Developmental_Trajectories_HCPD.pdf", device = "pdf", dpi = 500, width = 2.02 , height = 1.65)Region-specific GAM smooths
#Function for connection-specific GAM spline + derivative plots
connection_trajectory_plot <- function(measure, atlas, dataset, tract_name, display_color, y_ticks){
#Read in participant-level data and model-level fitted values and derivatives
participant.data.df <- get(sprintf("%s.%s.%s", measure, atlas, dataset))
fittedvalues.df <- get(sprintf("fittedvalues_%s_%s_%s", measure, atlas, dataset)) %>% filter(., tract == tract_name)
derivatives.df <- get(sprintf("derivatives_%s_%s_%s", measure, atlas, dataset)) %>% filter(., tract == tract_name)
derivatives.df$significant.derivative[derivatives.df$significant.derivative == 0] <- NA
#Connection spline plot with participant-level data
trajectory.plot <- ggplot(data = participant.data.df, aes(x = age, y = get(tract_name))) +
geom_point(color = alpha(display_color, 0.5), size = .1) +
geom_line(data = fittedvalues.df, aes(x = age, y = fitted), color = display_color, linewidth = 1) +
geom_ribbon(data = fittedvalues.df, aes(x = age, y = fitted, ymin = lower, ymax = upper), alpha = .8, linetype = 0, fill = display_color) +
labs(x="\nAge", y=sprintf("%s %s: %s\n", tract_name, measure, toupper(dataset))) +
theme_classic() +
scale_x_continuous(breaks = c(8, 12, 16, 20), limits = c(8,23), expand = c(0.025,.45)) +
scale_y_continuous(breaks = y_ticks) +
theme(
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2))
#Significant derivatives (rate of age-related change) plot
derivatives.plot <- ggplot(data = derivatives.df) +
geom_tile(aes(x = age, y = .1, fill = significant.derivative)) +
scale_fill_gradient(low = alpha(display_color, 0.2), high = display_color, na.value = "white") +
labs(x="\nAge", fill = "Derivative") +
theme_classic() +
scale_x_continuous(breaks=c(8, 12, 16, 20), limits = c(8,23), expand = c(0.025,.45)) +
theme(
legend.position = "none",
axis.title.y = element_blank(),
axis.text = element_blank(),
axis.line = element_blank(),
axis.ticks = element_blank(),
axis.title = element_blank(),
legend.text = element_text(size = 6, family = "Arial", color = c("black")))
allplots <- list(trajectory.plot, derivatives.plot)
tract.plot <- cowplot::plot_grid(rel_heights = c(16, 3.2), plotlist = allplots, align = "v", axis = "lr", ncol = 1)
return(tract.plot)
}Primary motor cortex (4)
Left 4: S-A axis rank = 16, mesulam assignment = idiotypic
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_4_autotrack", display_color = "#FEC22F", y_ticks = c(0.4, 0.45, 0.5, 0.55))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/PrimaryMotor_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_4_autotrack") %>% select(smooth.increase.offset)## smooth.increase.offset
## 1 16.06784
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_4_autotrack", display_color = "#FEC22F", y_ticks = c(0.35, 0.4, 0.45, 0.5))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/PrimaryMotor_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_4_autotrack") %>% select(smooth.increase.offset)## smooth.increase.offset
## 1 15.38233
Lateral parietal (PF)
Left PF: S-A axis rank = 248, mesulam assignment = heteromodal
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_PF_autotrack", display_color = "#672975", y_ticks = c(0.3, 0.35, 0.4, 0.45))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralParietal_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_PF_autotrack") %>% select(smooth.increase.offset)## smooth.increase.offset
## 1 18.4531
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_PF_autotrack", display_color = "#672975", y_ticks = c(0.3, 0.35, 0.4, 0.45))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralParietal_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_PF_autotrack") %>% select(smooth.increase.offset)## smooth.increase.offset
## 1 16.42504
Lateral prefrontal (8BL)
Left 8BL: S-A axis rank = 342, mesulam assignment = heteromodal
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "pnc", tract_name = "thalamus_L_8BL_autotrack", display_color = "#323280", y_ticks = c(0.3, 0.35, 0.4, 0.45))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralFrontal_PNC.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_pnc %>% filter(tract == "thalamus_L_8BL_autotrack") %>% select(smooth.increase.offset)## smooth.increase.offset
## 1 21.21106
connection_trajectory_plot(measure = "FA", atlas = "glasser", dataset = "hcpd", tract_name = "thalamus_L_8BL_autotrack", display_color = "#323280", y_ticks = c(0.3, 0.35, 0.4, 0.45))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/LateralFrontal_HCPD.pdf", device = "pdf", dpi = 500, width = 1.95, height = 1.8)
gameffects_FA_glasser_hcpd %>% filter(tract == "thalamus_L_8BL_autotrack") %>% select(smooth.increase.offset)## smooth.increase.offset
## 1 21.91667
Maturational patterns are highly reproducible
Correlation of thalamocortical connection development effect sizes between datasets
Spearman’s rho
gameffects.merged <- merge(gameffects_FA_glasser_pnc, gameffects_FA_glasser_hcpd, by = c("tract", "label", "orig_parcelname", "mesulam_assignment"), suffixes = c("_pnc", "_hcpd"))
cor.test(gameffects.merged$GAM.smooth.partialR2_pnc, gameffects.merged$GAM.smooth.partialR2_hcpd)##
## Pearson's product-moment correlation
##
## data: gameffects.merged$GAM.smooth.partialR2_pnc and gameffects.merged$GAM.smooth.partialR2_hcpd
## t = 15.505, df = 219, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.6538415 0.7808211
## sample estimates:
## cor
## 0.7233925
Spatial permutation (spin) based p-value
spin.data <- left_join(spin.parcels, gameffects.merged, by = c("label", "orig_parcelname", "tract")) #full set of parcel data in rh --> lh order for spinning. spin test null correlations use complete obs only. Each null correlation correspondence statistic is thus computed on a slightly reduced set of data, akin to a jackknife procedure
perm.sphere.p(x = spin.data$GAM.smooth.partialR2_pnc, y = spin.data$GAM.smooth.partialR2_hcpd, perm.id = perm.id.full, corr.type = "pearson") ## [1] 0
Correlation plot
ggplot(gameffects.merged, aes(x = GAM.smooth.partialR2_pnc, y = GAM.smooth.partialR2_hcpd)) +
geom_point(color = c("#FCAB6A"), size = 0.5) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
theme_classic() +
scale_x_continuous(limits = c(-0.01, 0.17)) +
scale_y_continuous(limits = c(-0.06, 0.15)) +
labs(x="\nPNC", y="HCPD\n") +
theme(
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Ageeffect_correlation_PNCHCPD.pdf", device = "pdf", dpi = 500, width = 1.99, height = 1.58)Correlation of thalamocortical connection age of maturation between datasets
Spearman’s rho
cor.test(gameffects.merged$smooth.increase.offset_pnc, gameffects.merged$smooth.increase.offset_hcpd)##
## Pearson's product-moment correlation
##
## data: gameffects.merged$smooth.increase.offset_pnc and gameffects.merged$smooth.increase.offset_hcpd
## t = 7.8687, df = 150, p-value = 6.581e-13
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.4172322 0.6442871
## sample estimates:
## cor
## 0.540529
Spatial permutation (spin) based p-value
perm.sphere.p(x = spin.data$smooth.increase.offset_pnc, y = spin.data$smooth.increase.offset_hcpd, perm.id = perm.id.full, corr.type = "pearson") ## [1] 0.00025
Correlation plot
ggplot(gameffects.merged, aes(x = smooth.increase.offset_pnc, y = smooth.increase.offset_hcpd)) +
geom_point(color = c("#9c3a80"), size = 0.4) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
theme_classic() +
scale_x_continuous(limits = c(14, 23), breaks = c(14, 16, 18, 20, 22), expand = c(0.05, 0)) +
scale_y_continuous(limits = c(14, 23), breaks = c(14, 16, 18, 20, 22), expand = c(0, 1)) +
labs(x="\nPNC", y="HCPD\n") +
theme(
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) Correlation of developmental magnitude and age of maturation within dataset
PNC
Spearman’s rho
cor.test(gameffects_FA_glasser_pnc$GAM.smooth.partialR2, gameffects_FA_glasser_pnc$smooth.increase.offset, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: gameffects_FA_glasser_pnc$GAM.smooth.partialR2 and gameffects_FA_glasser_pnc$smooth.increase.offset
## S = 229591, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.8278022
Spatial permutation (spin) based p-value
spin.data <- left_join(spin.parcels, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract"))
perm.sphere.p(x = spin.data$GAM.smooth.partialR2, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") ## [1] 0
Correlation plot
ggplot(gameffects_FA_glasser_pnc, aes(x = GAM.smooth.partialR2, y = smooth.increase.offset)) +
geom_point(color = c("#FCAB6A"), size = 0.5) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
theme_classic() +
scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
labs(x="\nAge effect (partial R2)", y="Age of maturation (years)\n") +
theme(
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/FA_Ageeffect_Ageofmaturation_PNC.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.58)HCPD
cor.test(gameffects_FA_glasser_hcpd$GAM.smooth.partialR2, gameffects_FA_glasser_hcpd$smooth.increase.offset, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: gameffects_FA_glasser_hcpd$GAM.smooth.partialR2 and gameffects_FA_glasser_hcpd$smooth.increase.offset
## S = 122759, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.8267486
Spatial permutation (spin) based p-value
spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract"))
perm.sphere.p(x = spin.data$GAM.smooth.partialR2, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") ## [1] 0
Early and Late Maturing Thalamocortical Connections
Early and late maturing areas
#Identify parcels that fall within first and fourth quartiles of the age of maturation variable
matstage.df <- gameffects_FA_glasser_pnc %>% select(label, smooth.increase.offset)
matstage.df$maturation_bracket <- "mid" #set default bracket to mid
matstage.df$maturation_bracket[matstage.df$smooth.increase.offset < quantile(matstage.df$smooth.increase.offset, 0.25, na.rm = T)] <- "early" #set first quartile bracket to early
matstage.df$maturation_bracket[matstage.df$smooth.increase.offset > quantile(matstage.df$smooth.increase.offset, 0.75, na.rm = T)] <- "late" #set fourth quartile bracket to late
matstage.df <- rbind(matstage.df, data.frame(label = "rh_???", smooth.increase.offset = NA, maturation_bracket = "medialwall")) #create a label for medial wall to color it grey
matstage.df <- rbind(matstage.df, data.frame(label = "lh_???", smooth.increase.offset = NA, maturation_bracket = "medialwall")) #create a label for medial wall to color it grey
ggseg(.data = matstage.df, atlas = "glasser", mapping = aes(fill = maturation_bracket, colour=I("black"), size=I(.06))) + scale_fill_manual(values = c("#FEC22F", alpha("#323280", 0.97), "white", "white" ), na.value = "gray90") + theme(legend.position = "none")## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi side region label smooth.increase.offset maturation_bracket
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L… 16.9 mid
NeuroSynth decoding of age of maturation maps
Compute cognitive term-development map spatial correlations
#Function to calculate the correlation between a neurosynth term z-score map and a thalamocortical development map
developmentmap.neurosynth.decoding <- function(df, dev.metric, term){
df.allregions <- left_join(spin.parcels, df, by = c("label", "tract", "orig_parcelname")) #full set of parcel data in rh --> lh order
term.cor <- cor(df.allregions[,term], df.allregions[,dev.metric], use = "complete.obs", method = c("spearman")) #correlation between neurosynth term map and development metric
term.results <- data.frame("term" = term, "term.correlation" = term.cor)
return(term.results)
}Developmental timing decoding by Neurosynth cognitive terms
PNC
gameffects_neurosynth_pnc <- merge(gameffects_FA_glasser_pnc, neurosynth.terms, by = "label")
devpattern.neurosynth.pnc <- ldply(neurosynth.termlist, function(t){
developmentmap.neurosynth.decoding(df = gameffects_neurosynth_pnc, dev.metric = "smooth.increase.offset", term = t)}) %>% #neurosynth-neurodevelopment correlations for all terms
arrange(desc(term.correlation))devpattern.neurosynth.pnc <- rbind(slice_max(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.pnc, order_by = term.correlation, n = 10)) #10 most positively correlated and negatively correlated neurosynth terms with the age of maturation map
devpattern.neurosynth.pnc## term term.correlation
## 1 expectancy 0.4111486
## 2 monitoring 0.3976832
## 3 cognitive_control 0.3815389
## 4 reasoning 0.3296662
## 5 strategy 0.3109975
## 6 response_inhibition 0.3094288
## 7 memory_retrieval 0.2866834
## 8 thought 0.2835596
## 9 memory 0.2775858
## 10 recall 0.2677131
## 11 object_recognition -0.2935350
## 12 visual_perception -0.2879558
## 13 morphology -0.2742436
## 14 consolidation -0.2560287
## 15 gaze -0.2521194
## 16 coordination -0.2340298
## 17 expertise -0.2327537
## 18 perception -0.2242119
## 19 facial_expression -0.2200723
## 20 induction -0.2038560
HCPD
gameffects_neurosynth_hcpd <- merge(gameffects_FA_glasser_hcpd, neurosynth.terms, by = "label")
devpattern.neurosynth.hcpd <- ldply(neurosynth.termlist, function(t){
developmentmap.neurosynth.decoding(df = gameffects_neurosynth_hcpd, dev.metric = "smooth.increase.offset", term = t)}) %>%
arrange(desc(term.correlation))devpattern.neurosynth.hcpd <- rbind(slice_max(devpattern.neurosynth.hcpd, order_by = term.correlation, n = 10), slice_min(devpattern.neurosynth.hcpd, order_by = term.correlation, n = 10))
devpattern.neurosynth.hcpd## term term.correlation
## 1 expectancy 0.3940717
## 2 sentence_comprehension 0.3852283
## 3 thought 0.3694649
## 4 cognitive_control 0.3561159
## 5 response_inhibition 0.3480972
## 6 monitoring 0.3193377
## 7 emotion_regulation 0.3186345
## 8 recall 0.2951983
## 9 language 0.2831427
## 10 inhibition 0.2653346
## 11 object_recognition -0.3456587
## 12 multisensory -0.3419028
## 13 visual_perception -0.3341608
## 14 adaptation -0.2941774
## 15 mental_imagery -0.2867927
## 16 induction -0.2677590
## 17 localization -0.2631644
## 18 fixation -0.2611707
## 19 gaze -0.2601696
## 20 imagery -0.2521325
Overlapping terms across datasets
merge(devpattern.neurosynth.pnc, devpattern.neurosynth.hcpd, by = "term") %>% arrange(desc(term.correlation.x)) %>% select(term)## term
## 1 expectancy
## 2 monitoring
## 3 cognitive_control
## 4 response_inhibition
## 5 thought
## 6 recall
## 7 induction
## 8 gaze
## 9 visual_perception
## 10 object_recognition
ggplot(devpattern.neurosynth.pnc, aes(x = term.correlation, y = term, color = term.correlation)) +
geom_segment(aes(y = reorder(term, term.correlation), yend = term, x = 0, xend =
term.correlation), color = "#989898", linewidth = 0.2) +
geom_point(size = 2.2, alpha = 0.85) +
theme_light() +
scale_color_gradientn(colors = c("#FEC22F", "#F59A72", "#9c3a80", "#672975","#323280"), limits = c(-0.29, 0.35), oob=squish) +
labs(x = "\nCognitive Term-Thalamocortical Development Correlation", y = "NeuroSynth Term\n") +
scale_x_continuous(limits = c(-0.3,0.42), breaks = c(-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4)) +
scale_y_discrete(
labels = c("expectancy"=expression(bold(expectancy)),
"monitoring"=expression(bold(monitoring)),
"cognitive_control"=expression(bold(cognitive_control)),
"response_inhibition"=expression(bold(response_inhibition)),
"thought"=expression(bold(thought)),
"recall"=expression(bold(recall)),
"induction"=expression(bold(induction)),
"object_recognition"=expression(bold(object_recognition)),
"visual_perception"=expression(bold(visual_perception)),
"gaze"=expression(bold(gaze)),
"induction"=expression(bold(induction), parse=TRUE))) +
theme(
legend.position = "none",
panel.grid.major.y = element_blank(),
panel.grid.major.x = element_line(color = c("gray90")),
axis.ticks.y = element_blank(),
panel.border = element_blank(),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.ticks = element_line(linewidth = .2))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure4/Neurosynth_Maturationdecoding.pdf", device = "pdf", dpi = 500, width = 3.35, height = 3.07)Probability of obtaining 10/20 overlapping words given the set of 123 terms
#Function to compute the number of overlapping terms when drawing term sets of size 20 at random for PNC and HCPD nulls
compute_term_overlap <- function(){
term.subset.pnc <- sample(x = neurosynth.termlist, size = 20, replace = FALSE)
term.subset.hcpd <- sample(x = neurosynth.termlist, size = 20, replace = FALSE)
num.overlap <- length(intersect(term.subset.pnc, term.subset.hcpd))
return(num.overlap)}
#Null distribution of overlapping term counts based on 10,000 random draws
set.seed(1)
overlap.number <- replicate(10000, compute_term_overlap()) %>% as.data.frame %>% set_names("overlaps")
sprintf("In %s out of 10,000 runs, the term overlap number was equal to or greater than 10", sum(overlap.number$overlaps >= 10))## [1] "In 1 out of 10,000 runs, the term overlap number was equal to or greater than 10"
sprintf("The permutation-based p value for these overlap counts is %s", sum(overlap.number$overlaps >= 10) / 10000)## [1] "The permutation-based p value for these overlap counts is 1e-04"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 2.000 3.000 3.257 4.000 10.000
ggplot(overlap.number, aes(x = overlaps)) +
geom_histogram(bins = 12, fill = c( "#9c3a80"), color = c("white")) +
theme_classic() +
scale_x_continuous(limits = c(-1, 10), breaks = c(0, 2, 4, 6, 8, 10)) +
labs(x = "\nNeurosynth term overlap", y = "Count\n") +
geom_vline(xintercept = 10, linewidth = .1) +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .1),
axis.ticks = element_line(linewidth = .1))Thalamocortical Connectivity Maturation is Organized by the Sensorimotor-Association Axis
gameffects_FA_glasser_pnc <- merge(gameffects_FA_glasser_pnc, (S.A.axis %>% select(tract, SA.axis)), by = "tract")
gameffects_FA_glasser_hcpd <- merge(gameffects_FA_glasser_hcpd, (S.A.axis %>% select(tract, SA.axis)), by = "tract")Age-resolved developmental change plots
PNC
derivatives_FA_glasser_pnc$significant.derivative[derivatives_FA_glasser_pnc$significant.derivative <= 0] <- NA #replace non-significantly increasing derivatives with NA. (Note here 0 was == not significant)
ggplot() +
geom_line(data = derivatives_FA_glasser_pnc, aes(x = age, y = SA.axis, group = tract, alpha = significant.derivative, color = SA.axis, linewidth = significant.derivative)) +
scale_alpha_continuous(range = c(0, 1), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
scale_linewidth_continuous(range = c(0, 4.5), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
scale_color_gradient2(low = "#ffc933", mid = "#e9dcf2", high = "#6f1282", guide = "colorbar", midpoint = 180) +
scale_x_continuous(breaks = c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .25)) +
scale_y_continuous(expand = c(0.03, 0)) +
theme_classic() +
ylab("Sensorimotor-Association Axis Rank\n") +
xlab("\nAge") +
theme(
legend.position = "none",
axis.line = element_line(size = .2),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.ticks = element_line(linewidth = .2))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Derivativeplot_pnc.pdf", device = "pdf", dpi = 500, width = 2.3, height = 3)HCPD
derivatives_FA_glasser_hcpd$significant.derivative[derivatives_FA_glasser_hcpd$significant.derivative <= 0] <- NA #replace non-significantly increasing derivatives with NA. (Note here 0 was == not significant)
ggplot() +
geom_line(data = derivatives_FA_glasser_hcpd, aes(x = age, y = SA.axis, group = tract, alpha = significant.derivative, color = SA.axis, linewidth = significant.derivative)) +
scale_alpha_continuous(range = c(0, 1), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
scale_linewidth_continuous(range = c(0, 4.5), breaks = c(0, 0.001, 0.002, 0.003, 0.004, 0.005, 0.006), limits = c(0.0005, 0.006)) +
scale_color_gradient2(low = "#ffc933", mid = "#e9dcf2", high = "#6f1282", guide = "colorbar", midpoint = 180) +
scale_x_continuous(breaks = c(8, 10, 12, 14, 16, 18, 20, 22), expand = c(0, .25)) +
scale_y_continuous(expand = c(0.03, 0)) +
theme_classic() +
ylab("Sensorimotor-Association Axis Rank\n") +
xlab("\nAge") +
theme(
legend.position = "none",
axis.line = element_line(size = .2),
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.ticks = element_line(linewidth = .2))Thalamocortical connections mature at progressively older ages across the S-A axis
PNC
Spearman’s rho
cor.test(gameffects_FA_glasser_pnc$SA.axis, gameffects_FA_glasser_pnc$smooth.increase.offset, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: gameffects_FA_glasser_pnc$SA.axis and gameffects_FA_glasser_pnc$smooth.increase.offset
## S = 685321, p-value = 2.986e-13
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.4859962
Spatial permutation (spin) based p-value
spin.data <- left_join(spin.parcels, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract"))
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") ## [1] 8e-04
Correlation plot
ggplot(gameffects_FA_glasser_pnc, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
geom_point(size = 0.5) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_color_gradient2(low = "goldenrod1", mid = "#ede4f5", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
theme_classic() +
scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
labs(x="\nS-A axis", y="Age of maturation (years)\n") +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_agematuration_pnc.pdf", device = "pdf", dpi = 500, width = 1.93, height = 1.62)Brain map
Thalamocortical Age of Maturation
agemat.df <- gameffects_FA_glasser_pnc %>% select(label, smooth.increase.offset)
agemat.df <- left_join((spin.data %>% select(label)), agemat.df, by = "label")
agemat.df$cortex <- "cortex"
agemat.df <- rbind(agemat.df, data.frame(label = "rh_???", smooth.increase.offset = NA, cortex = "medialwall")) #create a label for medial wall to color it grey
agemat.df <- rbind(agemat.df, data.frame(label = "lh_???", smooth.increase.offset = NA, cortex = "medialwall")) #create a label for medial wall to color it grey
ggplot() +
geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = smooth.increase.offset, colour=I("black"), size=I(.05))) +
theme_void() +
scale_color_gradient2(low = "goldenrod1", mid = "#c4b0d6", high = "#6f1282", guide = "colorbar", aesthetics = "fill", midpoint = 17.5, na.value = "gray93", limits = c(16, 19), oob = squish) +
new_scale_fill() +
geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = cortex, colour=I("black"), size=I(.05))) +
scale_fill_manual(values = c(alpha("white", 0), "white")) +
theme(
legend.text = element_text(size = 5, family = "Arial", color = c("black")),
legend.key.size = unit(1, 'mm'),
legend.title = element_text(size = 5, family = "Arial", color = c("black")))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Ageofmaturation_corticalplot_pnc.pdf", device = "pdf", dpi = 500, width = 5, height = 6.5) Sensorimotor-Association Axis
S.A.axis.plot <- S.A.axis[S.A.axis$tract %in% tractlist.pnc$tract,] %>% select(label, SA.axis)
ggplot() +
geom_brain(data = S.A.axis.plot, atlas = glasser, mapping = aes(fill = SA.axis, colour=I("black"), size=I(.05))) +
theme_void() +
scale_color_gradient2(low = "goldenrod1", mid = "#c4b0d6", high = "#6f1282", guide = "colorbar", aesthetics = "fill", midpoint = 180, na.value = "gray93") +
new_scale_fill() +
geom_brain(data = agemat.df, atlas = glasser, mapping = aes(fill = cortex, colour=I("black"), size=I(.05))) +
scale_fill_manual(values = c(alpha("white", 0), "white")) +
theme(
legend.text = element_text(size = 5, family = "Arial", color = c("black")),
legend.key.size = unit(1, 'mm'),
legend.title = element_text(size = 5, family = "Arial", color = c("black")))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAaxis_corticalplot_pnc.pdf", device = "pdf", dpi = 500, width = 4.7, height = 6.2) HCPD
Spearman’s rho
cor.test(gameffects_FA_glasser_hcpd$SA.axis, gameffects_FA_glasser_hcpd$smooth.increase.offset, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: gameffects_FA_glasser_hcpd$SA.axis and gameffects_FA_glasser_hcpd$smooth.increase.offset
## S = 342656, p-value = 2.008e-12
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.5164053
Spatial permutation (spin) based p-value
spin.data <- left_join(spin.parcels, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract"))
perm.sphere.p(x = spin.data$SA.axis, y = spin.data$smooth.increase.offset, perm.id = perm.id.full, corr.type = "spearman") ## [1] 0.00185
Correlation plot
ggplot(gameffects_FA_glasser_hcpd, aes(x = SA.axis, y = smooth.increase.offset, color = SA.axis)) +
geom_point(size = 0.5) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_color_gradient2(low = "goldenrod1", mid = "#ece4f2", high = "#6f1282", guide = "colorbar", aesthetics = "color", midpoint = 180) +
theme_classic() +
scale_y_continuous(limits = c(10.7, 22), breaks = c(12, 14, 16, 18, 20, 22)) +
labs(x="\nS-A axis", y="Age of maturation (years)\n") +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) S-A axis age resolved analysis
PNC
#Calculate the median derivative - S-A axis correlation value and its 95% credible interval at 100 ages
colnames(temporalSAcorr_FA_glasser_pnc) <- sprintf("draw%s", seq(from = 1, to = ncol(temporalSAcorr_FA_glasser_pnc))) #set column names; 100 age x 10,000 posterior draw correlations df
deriv.SAaxis.mediancorr <- apply(temporalSAcorr_FA_glasser_pnc, 1, function(x){median(x)}) #median correlation value at each of the 100 ages
cor.credible.intervals <- apply(temporalSAcorr_FA_glasser_pnc, 1, function(x){quantile(x, probs = c(0.025, 0.975))}) #compute the credible interval for the correlation value at each age based on all draws
cor.credible.intervals <- t(cor.credible.intervals)
cor.credible.intervals <- as.data.frame(cor.credible.intervals)
cor.credible.intervals$age <- as.numeric(seq(8, 23, length.out = 100))
cor.credible.intervals$SA.correlation <- deriv.SAaxis.mediancorr
colnames(cor.credible.intervals) <- c("lower","upper","age","SA.correlation")ggplot(cor.credible.intervals, aes(x = age, y = SA.correlation, ymin = lower, ymax = upper)) +
geom_ribbon(fill = c(alpha("#323280", 0.4))) +
geom_line(linewidth = .2, color = "black") +
labs(x="\nAge", y="Developmental Alignment to the S-A Axis\n") +
theme_classic() +
scale_x_continuous(breaks=c(8, 12, 16, 20), expand = c(0, .1)) +
scale_y_continuous(breaks=c(0, 0.1, 0.2, 0.3, 0.4), limits = c(0, 0.4)) +
theme(
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Derivative-SAcorrelation_plot_pnc.pdf", device = "pdf", dpi = 500, width = 1.35, height = 1.03)HCPD
#Calculate the median derivative - S-A axis correlation value and its 95% credible interval at 100 ages
colnames(temporalSAcorr_FA_glasser_hcpd) <- sprintf("draw%s", seq(from = 1, to = ncol(temporalSAcorr_FA_glasser_hcpd))) #set column names; 100 age x 10,000 posterior draw correlations df
deriv.SAaxis.mediancorr <- apply(temporalSAcorr_FA_glasser_hcpd, 1, function(x){median(x)}) #median correlation value at each of the 100 ages
cor.credible.intervals <- apply(temporalSAcorr_FA_glasser_hcpd, 1, function(x){quantile(x, probs = c(0.025, 0.975))}) #compute the credible interval for the correlation value at each age based on all draws
cor.credible.intervals <- t(cor.credible.intervals)
cor.credible.intervals <- as.data.frame(cor.credible.intervals)
cor.credible.intervals$age <- as.numeric(seq(8, 23, length.out = 100))
cor.credible.intervals$SA.correlation <- deriv.SAaxis.mediancorr
colnames(cor.credible.intervals) <- c("lower","upper","age","SA.correlation")ggplot(cor.credible.intervals, aes(x = age, y = SA.correlation, ymin = lower, ymax = upper)) +
geom_ribbon(fill = c(alpha("#323280", 0.4))) +
geom_line(linewidth = .2, color = "black") +
labs(x="\nAge", y="Developmental Alignment to the S-A Axis\n") +
theme_classic() +
scale_x_continuous(breaks=c(8, 12, 16, 20), expand = c(0, .1)) +
#scale_y_continuous(breaks=c(0, 0.1, 0.2, 0.3, 0.4), limits = c(0, 0.4)) +
theme(
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) Thalamocortical maturation along anatomical axes
#Glasser parcel x,y,z coordinates
load(file = "/cbica/projects/thalamocortical_development/Maps/parcellations/surface/hcp_mmp1.0.rda") #RDA file from the brainGraph package https://github.com/cwatson/brainGraph/blob/master/data/hcp_mmp1.0.rda with glasser region x, y, and z MNI coordinates. regions are in lh --> rh order
hcp_mmp1.0$label <- NA
hcp_mmp1.0$label[1:180] <- glasser.regions$label[181:360]
hcp_mmp1.0$label[181:360] <- glasser.regions$label[1:180]#Function to compute the correlation between a development map and an anatomical axis (x,y,z) as well as test the significance of the difference between the developmental correlation with the anatomical axis versus the S-A axis
compute_axis_correlation <- function(df.dev, dev.metric, axis){
df.dev.axes <- merge(df.dev, hcp_mmp1.0, by = "label")
df.dev.axes.spin <- left_join(spin.parcels, df.dev.axes, by = c("label", "tract", "orig_parcelname"))
#S-A axis - development correlation (as computed above) for comparison of two overlapping correlations based on dependent groups
SA.axis.cor <- cor(df.dev.axes$SA.axis, df.dev.axes[,dev.metric], use = "complete.obs", method = c("spearman"))
#Anatomical axis - development correlation
anatomical.axis.cor <- cor(df.dev.axes[,axis], df.dev.axes[,dev.metric], use = "complete.obs", method = c("spearman")) #correlation between anatomical axis and development metric
anatomical.axis.pvalue <- perm.sphere.p(x = df.dev.axes.spin[,axis], y = df.dev.axes.spin[,dev.metric], perm.id = perm.id.full, corr.type = "spearman") #spin based p-value for the correlation
#S-A axis - anatomical axis correlation
SA.anatomical.cor <- cor(df.dev.axes$SA.axis, df.dev.axes[,axis], use = "complete.obs", method = c("spearman"))
#Test for significant difference between two correlations based on dependent groups (i.e., correlations with an overlapping input)
cocor.result <- cocor.pvalue <- cocor.dep.groups.overlap(
r.jk = SA.axis.cor,
r.jh = anatomical.axis.cor,
r.kh = SA.anatomical.cor,
n = nrow(df.dev.axes), alternative = "two.sided", test = "hittner2003")
cocor.pvalue <- cocor.result@hittner2003[["p.value"]]
comparison.results <- data.frame(axis = axis, axis.cor = anatomical.axis.cor, axis.cor.pvalue = anatomical.axis.pvalue, SA.cor = SA.axis.cor, cocor.pvalue = cocor.pvalue)
return(comparison.results)
}PNC
Anterior-posterior axis
compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "y.mni")## axis axis.cor axis.cor.pvalue SA.cor cocor.pvalue
## 1 y.mni 0.3097694 0.06185 0.4859962 0.0001762558
Medial-lateral axis
compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "x.mni")## axis axis.cor axis.cor.pvalue SA.cor cocor.pvalue
## 1 x.mni -0.09059357 0.56025 0.4859962 1.604827e-10
Dorsal-ventral axis
compute_axis_correlation(df.dev = gameffects_FA_glasser_pnc, dev.metric = "smooth.increase.offset", axis = "z.mni")## axis axis.cor axis.cor.pvalue SA.cor cocor.pvalue
## 1 z.mni 0.1192805 0.2582 0.4859962 0.000236177
FDR-corrected p-values across cocor comparisons
anatomical.ps <- c(0.0001762558, 1.604827e-10, 0.000236177)
p.adjust(anatomical.ps, method = c("fdr"))## [1] 2.361770e-04 4.814481e-10 2.361770e-04
#Put above results into df for plotting
axis.results <- data.frame(Axis = factor(), corr = double())
axis.results <- axis.results %>% add_row(Axis = "S-A", corr = 0.4859962)
axis.results <- axis.results %>% add_row(Axis = "A-P", corr = 0.3097694)
axis.results <- axis.results %>% add_row(Axis = "D-V", corr = 0.1192805)
axis.results <- axis.results %>% add_row(Axis = "M-L", corr = -0.09059357)
axis.results$Axis <- factor(axis.results$Axis, ordered = TRUE, levels = c("S-A", "A-P", "D-V", "M-L"))
ggplot(axis.results, aes(x = Axis, y = corr, fill = Axis)) +
geom_bar(stat='identity') +
labs(x="") +
labs(y="\nAge of Maturation Correlation") +
theme_classic()+
scale_fill_manual(values=c(alpha("#323280", 0.5), "#672975", "#672975", "#672975")) +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2))ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/Anatomicalaxes_agematuration_comparison_pnc.pdf", device = "pdf", dpi = 500, width = 1.4, height = 1)HCPD
Anterior-posterior axis
compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "y.mni")## axis axis.cor axis.cor.pvalue SA.cor cocor.pvalue
## 1 y.mni 0.47926 0.01445 0.5164053 0.3912361
Medial-lateral axis
compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "x.mni")## axis axis.cor axis.cor.pvalue SA.cor cocor.pvalue
## 1 x.mni -0.04929832 0.72215 0.5164053 7.947443e-11
Dorsal-ventral axis
compute_axis_correlation(df.dev = gameffects_FA_glasser_hcpd, dev.metric = "smooth.increase.offset", axis = "z.mni")## axis axis.cor axis.cor.pvalue SA.cor cocor.pvalue
## 1 z.mni 0.06423155 0.4144 0.5164053 4.340955e-06
FDR-corrected p-values across cocor comparisons
anatomical.ps <- c(0.3912361, 7.947443e-11, 4.340955e-06)
p.adjust(anatomical.ps, method = c("fdr"))## [1] 3.912361e-01 2.384233e-10 6.511432e-06
Thalamocortical Connectivity Variation Developmentally Increases in Association Cortex
#Function to compute the coefficient of variation across thalamocortical connections within bins of the S-A axis
compute_SAbin_cov <- function(df.dwi, tractlist){
df.dwi.long <- df.dwi %>% melt(., id.vars=c("rbcid", "age", "sex", "mean_fd"), variable.name = "tract", value.name = "dwi.metric") %>% merge(., S.A.axis, by = "tract", sort = F)
df.dwi.long <- df.dwi.long[df.dwi.long$tract %in% tractlist$tract,] #only consider tracts in this dataset's final tractlist
cov.SAbin <- df.dwi.long %>% group_by(rbcid, SA.axis.bin) %>% #for every individual and every bin of the S-A axis
do(cov = (cv(as.numeric(.$dwi.metric), na.rm = T))) %>% #calculate the coefficient of variation for the dwi metric within each bin
unnest(cols = c(cov))
cov.SAbin <- merge(cov.SAbin, (df.dwi %>% select(rbcid, age, sex, mean_fd)), by = "rbcid")
return(cov.SAbin) #df with rows length of dataset N*5 and rows of rbcid, SA.axis.bin, cov
}PNC
#gam interaction test to determine whether changes in cov by age vary across S-A axis bins
summary(gam(cov ~ s(age, k = 4) + s(SA.axis.bin, k = 4) + ti(age, SA.axis.bin, k = 4), method = "REML", data = cov.SAbins.pnc))##
## Family: gaussian
## Link function: identity
##
## Formula:
## cov ~ s(age, k = 4) + s(SA.axis.bin, k = 4) + ti(age, SA.axis.bin,
## k = 4)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1315122 0.0002498 526.4 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(age) 1.001 1.002 12.68 0.000371 ***
## s(SA.axis.bin) 2.998 3.000 3191.64 < 2e-16 ***
## ti(age,SA.axis.bin) 3.116 3.585 13.74 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.627 Deviance explained = 62.8%
## -REML = -14561 Scale est. = 0.00035727 n = 5725
gamresults.bins <- vector("list", length = 5)
for(bin in 1:5){
cov.SAbin.pnc <- cov.SAbins.pnc %>% filter(SA.axis.bin == bin)
gamresults.bins[[bin]] <- gam.fit.smooth(measure = "cov", atlas = "SAbin", dataset = "pnc", region = "cov", smooth_var = "age", covariates = "sex + mean_fd", knots = 3, set_fx = TRUE) %>% select(-tract) %>% mutate(SA.bin = bin)
}
gamresults.cov.pnc <- do.call(rbind, gamresults.bins)
gamresults.cov.pnc <- gamresults.cov.pnc %>% mutate(fdr.pvalue = p.adjust(Anova.smooth.pvalue, method = c("fdr")))
knitr::kable(gamresults.cov.pnc)| GAM.smooth.Fvalue | GAM.smooth.pvalue | GAM.smooth.partialR2 | Anova.smooth.pvalue | smooth.change.onset | smooth.peak.change | smooth.decrease.onset | smooth.decrease.offset | smooth.increase.onset | smooth.increase.offset | smooth.last.change | SA.bin | fdr.pvalue |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 8.8422351 | 0.0001546 | -0.0152757 | 0.0001546 | 8.166667 | 8.576633 | 8.166667 | 16.06784 | NA | NA | 15.99330 | 1 | 0.0003866 |
| 0.0355421 | 0.9650831 | 0.0000624 | 0.9650831 | NA | NA | NA | NA | NA | NA | NA | 2 | 0.9650831 |
| 8.4049572 | 0.0002379 | 0.0145313 | 0.0002379 | 11.371859 | 11.483668 | NA | NA | 11.371859 | 16.51508 | 16.44054 | 3 | 0.0003965 |
| 60.9790643 | 0.0000000 | 0.0966420 | 0.0000000 | 8.166667 | 8.762982 | NA | NA | 8.166667 | 19.19849 | 19.12395 | 4 | 0.0000000 |
| 3.5096789 | 0.0302301 | 0.0061197 | 0.0302301 | 14.651591 | 15.844221 | NA | NA | 14.651591 | 15.91876 | 15.84422 | 5 | 0.0377876 |
#Function to plot participant-level cov by age in a specific bin of the S-A axis
plot_SAbin_cov <- function(SA.bin, plotcolor){
cov.SAbin.pnc <- cov.SAbins.pnc %>% filter(SA.axis.bin == SA.bin)
cov_plot <- ggplot(cov.SAbin.pnc, aes(x = age, y = cov)) +
geom_point(color = alpha(plotcolor, 0.4), size = 0.5) +
geom_smooth(method = "gam", formula = y ~ s(x, k = 3, fx = T), linewidth = 0.2, color = "black") +
theme_classic() +
labs(x="\nAge", y="Across-tract variation in FA (cov)\n") +
scale_x_continuous(breaks = c(8, 12, 16, 20), expand = c(0,.45)) +
scale_y_continuous(breaks = c(0.04, 0.08, 0.12, 0.16, 0.20, 0.24)) +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2))
return(cov_plot)
}ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAbin1_cov_development_pnc.pdf", device = "pdf", dpi = 500, width = 1.67, height = 1.3)ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAbin2_cov_development_pnc.pdf", device = "pdf", dpi = 500, width = 1.67, height = 1.3)ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAbin3_cov_development_pnc.pdf", device = "pdf", dpi = 500, width = 1.67, height = 1.3)ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAbin4_cov_development_pnc.pdf", device = "pdf", dpi = 500, width = 1.67, height = 1.3)ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure5/SAbin5_cov_development_pnc.pdf", device = "pdf", dpi = 500, width = 1.67, height = 1.3)HCPD
#gam interaction test to determine whether changes in cov by age vary across S-A axis bins
summary(gam(cov ~ s(age, k = 4) + s(SA.axis.bin, k = 4) + ti(age, SA.axis.bin), method = "REML", data = cov.SAbins.hcpd))##
## Family: gaussian
## Link function: identity
##
## Formula:
## cov ~ s(age, k = 4) + s(SA.axis.bin, k = 4) + ti(age, SA.axis.bin)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1139555 0.0002661 428.3 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(age) 1.404 1.685 6.056 0.00296 **
## s(SA.axis.bin) 2.996 3.000 1895.366 < 2e-16 ***
## ti(age,SA.axis.bin) 6.540 8.427 6.289 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.668 Deviance explained = 66.9%
## -REML = -8065.3 Scale est. = 0.0002025 n = 2860
gamresults.bins <- vector("list", length = 5)
for(bin in 1:5){
cov.SAbin.hcpd <- cov.SAbins.hcpd %>% filter(SA.axis.bin == bin)
gamresults.bins[[bin]] <- gam.fit.smooth(measure = "cov", atlas = "SAbin", dataset = "hcpd", region = "cov", smooth_var = "age", covariates = "sex + mean_fd", knots = 3, set_fx = TRUE) %>% select(-tract) %>% mutate(SA.bin = bin)
}
gamresults.cov.hcpd <- do.call(rbind, gamresults.bins)
gamresults.cov.hcpd <- gamresults.cov.hcpd %>% mutate(fdr.pvalue = p.adjust(Anova.smooth.pvalue, method = c("fdr")))
knitr::kable(gamresults.cov.hcpd)| GAM.smooth.Fvalue | GAM.smooth.pvalue | GAM.smooth.partialR2 | Anova.smooth.pvalue | smooth.change.onset | smooth.peak.change | smooth.decrease.onset | smooth.decrease.offset | smooth.increase.onset | smooth.increase.offset | smooth.last.change | SA.bin | fdr.pvalue |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0.762425 | 0.4670114 | -0.0026821 | 0.4670114 | NA | NA | NA | NA | NA | NA | NA | 1 | 0.4670114 |
| 0.858181 | 0.4244824 | -0.0030180 | 0.4244824 | NA | NA | NA | NA | NA | NA | NA | 2 | 0.4670114 |
| 24.451553 | 0.0000000 | 0.0794006 | 0.0000000 | 8.083333 | 8.535176 | NA | NA | 8.083333 | 17.12018 | 17.05067 | 3 | 0.0000000 |
| 28.500997 | 0.0000000 | 0.0913491 | 0.0000000 | 8.083333 | 8.569933 | NA | NA | 8.083333 | 17.05067 | 16.98116 | 4 | 0.0000000 |
| 1.097592 | 0.3343815 | 0.0038566 | 0.3343815 | NA | NA | NA | NA | NA | NA | NA | 5 | 0.4670114 |
Thalamocortical Connectivity Maturation Mirrors Timescales of Cortical Plasticity
Plasticity marker development data
#Fluctuation amplitude age of decrease onset data from Sydnor et al., 2023, Nat Neurosci https://www.nature.com/articles/s41593-023-01282-y
boldamplitude.agedecrease.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/spatiotemp_dev_plasticity/cortical_maps/AgeofDeclineOnset_FirstNegDeriv.pscalar.nii")
boldamplitude.dev <- data.frame(orig_parcelname = names(boldamplitude.agedecrease.cifti$Parcel), amplitude.agedecline = boldamplitude.agedecrease.cifti$data)#T1/T2 ratio rate of developmental change from Baum et al., 2022, J Neuro https://www.jneurosci.org/content/42/29/5681
myelin.maxdev.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_median_posterior_age_of_max_slope_myelination.pscalar.nii") #age of maximal T1/T2 ratio increase
myelin.roc.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/myelin_development/hcpd_n628_mean_posterior_annualized_roc_myelin.pscalar.nii") #annualized rate of change
myelin.dev <- data.frame(orig_parcelname = names(myelin.maxdev.cifti$Parcel), myelin.agemaxdev = myelin.maxdev.cifti$data, myelin.annualroc = myelin.roc.cifti$data)#Generate the E:I ratio-age slope map in Glasser HCP-MMP parcels from Zhang, Larsen et al., bioRxiv, https://www.biorxiv.org/content/10.1101/2023.06.22.546023v1.full
library(ciftiTools)
#Calculate E:I age slopes in Yan100 parcels
EI.group.ages <- read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/age_Yan.csv", header = F) %>% as.data.frame %>% set_names("age") %>% mutate(age = age/12) #average age of 29 age groups
EI.group.ratio <- t(read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/EI_Yan.csv", header = F)) %>% as.data.frame() #average E:I ratio in all of the Yan100 parcels in the 29 age groups
Yan100.EI.ageslopes <- sapply(EI.group.ratio, function(x) lm(x ~ EI.group.ages$age)$coefficients[2]) %>% as.data.frame() %>% set_names("slope") #linear model coefficients for the association between age and E:I ratio in all Yan100 parcels
#Reorder Yan100 parcels to match the released version of the atlas
Yan100.region.reordering <- read.csv("/cbica/projects/thalamocortical_development/Maps/EI_development/Yanatlas_regionmapping.csv") #mapping of regions between version
Yan100.EI.ageslopes$mapping <- Yan100.region.reordering$parcelnumber.released
Yan100.EI.ageslopes <- Yan100.EI.ageslopes %>% arrange(mapping)
#Map parcel-level Yan100 data to fslr 32k vertex-level data
Yan100.densecifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/EI_development/100Parcels_Yeo2011_17Networks.dscalar.nii")
Yan100.lh.vertices <- Yan100.densecifti$data$cortex_left %>% as.data.frame %>% set_names("mapping") #parcel number assignments for 32492 vertices
Yan100.lh.vertices.slope <- left_join(Yan100.lh.vertices, Yan100.EI.ageslopes, by = "mapping") %>% select("slope") #age slope for all vertices
Yan100.rh.vertices <- Yan100.densecifti$data$cortex_right %>% as.data.frame %>% set_names("mapping") #parcel number assignments for 32492 vertices
Yan100.rh.vertices.slope <- left_join(Yan100.rh.vertices, Yan100.EI.ageslopes, by = "mapping") %>% select("slope") #age slope for all vertices
lh.medialwall <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/medialwall.mask.leftcortex.csv", header = F)
rh.medialwall <- read.csv("/cbica/projects/thalamocortical_development/Maps/parcellations/surface/medialwall.mask.rightcortex.csv", header = F)
Yan100.EI.ageslopes.densecifti <- as_cifti(cortexL = Yan100.lh.vertices.slope$slope, cortexR = Yan100.rh.vertices.slope$slope, cortexL_mwall = lh.medialwall$V1, cortexR_mwall = rh.medialwall$V1) #ciftify me captain
write_cifti(Yan100.EI.ageslopes.densecifti, "/cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_Yan100.dscalar.nii")
#Parcellate the fslr 32k vertex-level E:I-age slope with the glasser HCP-MMP atlas
command = "-cifti-parcellate /cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_Yan100.dscalar.nii /cbica/projects/thalamocortical_development/Maps/parcellations/surface/glasser_space-fsLR_den-32k_desc-atlas.dlabel.nii COLUMN /cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_glasser360.pscalar.nii -only-numeric"
ciftiTools::run_wb_cmd(command, intern = FALSE, ignore.stdout = NULL, ignore.stderr = NULL)
detach("package:ciftiTools", unload=TRUE)EI.ageslope.cifti <- read_cifti("/cbica/projects/thalamocortical_development/Maps/EI_development/EI_ageslopes_glasser360.pscalar.nii")
EIratio.dev <- data.frame(orig_parcelname = names(EI.ageslope.cifti$Parcel), EI.ageslope = EI.ageslope.cifti$data)df.list <- list(boldamplitude.dev, myelin.dev, EIratio.dev) #dfs to merge
plasticity.dev <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list)
plasticity.dev <- left_join(spin.parcels, plasticity.dev, by = "orig_parcelname")Fluctuation amplitude decrease onset map
ggplot() +
geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = amplitude.agedecline, colour=I("black"), size=I(.03))) +
theme_void() +
scale_fill_gradientn(colors = c("#fcd16d", "#EA8B57","#943b90", "#581768"), limits = c(7.9, 18), oob = squish, na.value = "white") ## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi region side label geometry orig_parcelname
## <chr> <chr> <chr> <chr> <chr> <chr> <MULTIPOLYGON> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L… EMPTY L_10pp_ROI
## # ℹ 5 more variables: tract <chr>, amplitude.agedecline <dbl>,
## # myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
## label orig_parcelname tract amplitude.agedecline
## 396 lh_L_10pp L_10pp_ROI thalamus_L_10pp_autotrack 17.18593
## myelin.agemaxdev myelin.annualroc EI.ageslope
## 396 16.37678 0.008886078 -0.03484702
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Boldamplitude_agedecline_corticalmap.pdf", device = "pdf", dpi = 500, width = 4.75, height = 5.4) ## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi region side label geometry orig_parcelname
## <chr> <chr> <chr> <chr> <chr> <chr> <MULTIPOLYGON> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L… EMPTY L_10pp_ROI
## # ℹ 5 more variables: tract <chr>, amplitude.agedecline <dbl>,
## # myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
## label orig_parcelname tract amplitude.agedecline
## 396 lh_L_10pp L_10pp_ROI thalamus_L_10pp_autotrack 17.18593
## myelin.agemaxdev myelin.annualroc EI.ageslope
## 396 16.37678 0.008886078 -0.03484702
T1/T2 annualized rate of increase map
ggplot() +
geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = myelin.annualroc, colour=I("black"), size=I(.03))) +
theme_void() +
scale_fill_gradientn(colors = c( "#581768", "#943b90", "#EA8B57", "#fcd16d"), limits = c(0.005, 0.018), oob = squish, na.value = "white") ## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi region side label geometry orig_parcelname
## <chr> <chr> <chr> <chr> <chr> <chr> <MULTIPOLYGON> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L… EMPTY L_10pp_ROI
## # ℹ 5 more variables: tract <chr>, amplitude.agedecline <dbl>,
## # myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
## label orig_parcelname tract amplitude.agedecline
## 396 lh_L_10pp L_10pp_ROI thalamus_L_10pp_autotrack 17.18593
## myelin.agemaxdev myelin.annualroc EI.ageslope
## 396 16.37678 0.008886078 -0.03484702
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2ratio_annualizedroc_corticalmap.pdf", device = "pdf", dpi = 500, width = 4.44, height = 5.13) ## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi region side label geometry orig_parcelname
## <chr> <chr> <chr> <chr> <chr> <chr> <MULTIPOLYGON> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L… EMPTY L_10pp_ROI
## # ℹ 5 more variables: tract <chr>, amplitude.agedecline <dbl>,
## # myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
## label orig_parcelname tract amplitude.agedecline
## 396 lh_L_10pp L_10pp_ROI thalamus_L_10pp_autotrack 17.18593
## myelin.agemaxdev myelin.annualroc EI.ageslope
## 396 16.37678 0.008886078 -0.03484702
E:I ratio magnitude of developmental decline map
ggplot() +
geom_brain(data = plasticity.dev, atlas = glasser, mapping = aes(fill = EI.ageslope, colour=I("black"), size=I(.03))) +
theme_void() +
scale_fill_gradientn(colors = c("#fcd16d", "#EA8B57","#943b90", "#581768"), limits = c(-0.04, -0.0345), oob = squish, na.value = "white") ## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi region side label geometry orig_parcelname
## <chr> <chr> <chr> <chr> <chr> <chr> <MULTIPOLYGON> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L… EMPTY L_10pp_ROI
## # ℹ 5 more variables: tract <chr>, amplitude.agedecline <dbl>,
## # myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
## label orig_parcelname tract amplitude.agedecline
## 396 lh_L_10pp L_10pp_ROI thalamus_L_10pp_autotrack 17.18593
## myelin.agemaxdev myelin.annualroc EI.ageslope
## 396 16.37678 0.008886078 -0.03484702
ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIratio_decline_corticalmap.pdf", device = "pdf", dpi = 500, width = 4.13, height = 5.08) ## merging atlas and data by 'label'
## Warning: Some data not merged properly. Check for naming errors in data:
## atlas type hemi region side label geometry orig_parcelname
## <chr> <chr> <chr> <chr> <chr> <chr> <MULTIPOLYGON> <chr>
## 1 <NA> <NA> <NA> <NA> <NA> lh_L… EMPTY L_10pp_ROI
## # ℹ 5 more variables: tract <chr>, amplitude.agedecline <dbl>,
## # myelin.agemaxdev <dbl>, myelin.annualroc <dbl>, EI.ageslope <dbl>
## Warning: Some data not merged. Check for spelling mistakes in:
## label orig_parcelname tract amplitude.agedecline
## 396 lh_L_10pp L_10pp_ROI thalamus_L_10pp_autotrack 17.18593
## myelin.agemaxdev myelin.annualroc EI.ageslope
## 396 16.37678 0.008886078 -0.03484702
Thalamocortical maturational patterns align to the spatiotemporal unfolding of plasticity readouts
PNC
plasticity.dev.pnc <- left_join(plasticity.dev, gameffects_FA_glasser_pnc, by = c("label", "orig_parcelname", "tract"))Thalamocortical connectivity maturation - E:I ratio magnitude of developmental decline
Spearman’s rho
cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$EI.ageslope, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: plasticity.dev.pnc$smooth.increase.offset and plasticity.dev.pnc$EI.ageslope
## S = 735977, p-value = 2.893e-11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.4480033
Spatial permutation (spin) based p-value
EI.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$EI.ageslope, perm.id = perm.id.full, corr.type = "spearman") Correlation plot
ggplot(plasticity.dev.pnc, aes(x = EI.ageslope, y = smooth.increase.offset)) +
geom_point(size = 0.5, color = "#fcd16d") +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
theme_classic() +
labs(x="E/I ratio decline (age slope)\n", y="\nThalamocortical age of maturation (years)") +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIdecline_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) Thalamocortical connectivity maturation - myelin annualized rate of growth
Spearman’s rho
cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$myelin.annualroc, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: plasticity.dev.pnc$smooth.increase.offset and plasticity.dev.pnc$myelin.annualroc
## S = 1927545, p-value = 3.753e-11
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.4456946
Spatial permutation (spin) based p-value
myelin.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$myelin.annualroc, perm.id = perm.id.full, corr.type = "spearman") Correlation plot
ggplot(plasticity.dev.pnc, aes(x = myelin.annualroc, y = smooth.increase.offset, color = myelin.annualroc)) +
geom_point(size = 0.5, color = "#F1A45F") +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_x_continuous(limits = c(0.005, 0.025)) +
scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
theme_classic() +
labs(x="T1/T2 ratio annualized growth rate\n", y="\nThalamocortical age of maturation (years)", ) +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2growth_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) Thalamocortical connectivity maturation - fluctuation amplitude age of decrease onset
Spearman’s rho
cor.test(plasticity.dev.pnc$smooth.increase.offset, plasticity.dev.pnc$amplitude.agedecline, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: plasticity.dev.pnc$smooth.increase.offset and plasticity.dev.pnc$amplitude.agedecline
## S = 710935, p-value = 2.877e-05
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.303949
Spatial permutation (spin) based p-value
amplitude.thalamus.pnc <- perm.sphere.p(x = plasticity.dev.pnc$smooth.increase.offset, y = plasticity.dev.pnc$amplitude.agedecline, perm.id = perm.id.full, corr.type = "spearman")Correlation plot
ggplot(plasticity.dev.pnc, aes(x = amplitude.agedecline, y = smooth.increase.offset, color = amplitude.agedecline)) +
geom_point(size = 0.5, color = "#7B2C80") +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_x_continuous(limits = c(8.2, 21.2), breaks = c(8, 10, 12, 14, 16, 18, 20)) +
scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
theme_classic() +
labs(x="Amplitude decrease onset (years)\n", y="\nThalamocortical age of maturation (years)") +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Amplitudedecline_thalamocortical_pnc.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) FDR-corrected p-values across correlations
## [1] 0.00025 0.00085 0.02850
## [1] 0.000750 0.001275 0.028500
HCPD
plasticity.dev.hcpd <- left_join(plasticity.dev, gameffects_FA_glasser_hcpd, by = c("label", "orig_parcelname", "tract"))Thalamocortical connectivity maturation - E:I ratio magnitude of developmental decline
Spearman’s rho
cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$EI.ageslope, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: plasticity.dev.hcpd$smooth.increase.offset and plasticity.dev.hcpd$EI.ageslope
## S = 388563, p-value = 1.62e-09
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.4516166
Spatial permutation (spin) based p-value
EI.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$EI.ageslope, perm.id = perm.id.full, corr.type = "spearman") Correlation plot
ggplot(plasticity.dev.hcpd, aes(x = EI.ageslope, y = smooth.increase.offset, color = EI.ageslope)) +
geom_point(size = 0.5, color = "#fcd16d") +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
theme_classic() +
labs(x="E/I ratio decline (age slope)\n", y="\nThalamocortical age of maturation (years)") +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/EIdecline_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) Thalamocortical connectivity maturation - myelin annualized rate of growth
Spearman’s rho
cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$myelin.annualroc, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: plasticity.dev.hcpd$smooth.increase.offset and plasticity.dev.hcpd$myelin.annualroc
## S = 1032614, p-value = 9.472e-10
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.4573402
Spatial permutation (spin) based p-value
myelin.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$myelin.annualroc, perm.id = perm.id.full, corr.type = "spearman") Correlation plot
ggplot(plasticity.dev.hcpd, aes(x = myelin.annualroc, y = smooth.increase.offset, color = myelin.annualroc)) +
geom_point(size = 0.5, color = "#F1A45F") +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_x_continuous(limits = c(0.005, 0.025)) +
scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
theme_classic() +
labs(x="T1/T2 ratio annualized growth rate\n", y="\nThalamocortical age of maturation (years)", ) +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/T1T2growth_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) Thalamocortical connectivity maturation - fluctuation amplitude age of decrease onset
Spearman’s rho
cor.test(plasticity.dev.hcpd$smooth.increase.offset, plasticity.dev.hcpd$amplitude.agedecline, method = c("spearman"), exact = F)##
## Spearman's rank correlation rho
##
## data: plasticity.dev.hcpd$smooth.increase.offset and plasticity.dev.hcpd$amplitude.agedecline
## S = 346503, p-value = 8.766e-09
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.4416828
Spatial permutation (spin) based p-value
amplitude.thalamus.hcpd <- perm.sphere.p(x = plasticity.dev.hcpd$smooth.increase.offset, y = plasticity.dev.hcpd$amplitude.agedecline, perm.id = perm.id.full, corr.type = "spearman") Correlation plot
ggplot(plasticity.dev.hcpd, aes(x = amplitude.agedecline, y = smooth.increase.offset, color = amplitude.agedecline)) +
geom_point(size = 0.5, color = "#7B2C80") +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
scale_x_continuous(limits = c(8.2, 21.2), breaks = c(8, 10, 12, 14, 16, 18, 20)) +
scale_y_continuous(limits = c(12, 23), breaks = c(12, 14, 16, 18, 20, 22)) +
theme_classic() +
labs(x="Amplitude decrease onset (years)\n", y="\nThalamocortical age of maturation (years)") +
theme(
legend.position = "none",
axis.text = element_text(size = 5, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 5, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 5, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) ggsave(filename = "/cbica/projects/thalamocortical_development/figures/images/Figure6/Amplitudedecline_thalamocortical_hcpd.pdf", device = "pdf", dpi = 500, width = 1.8, height = 1.35) FDR-corrected p-values across correlations
## [1] 0.00020 0.00115 0.00590
## [1] 0.000600 0.001725 0.005900